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Abstract 


The  Cosmic  Ray  Inspection  and  Passive  Tomography  (CRIPT)  collaboration  has  con¬ 
structed  a  large-scale  detector  prototype  for  investigating  the  use  of  comic  ray  muon  scat¬ 
tering  tomography  for  Special  Nuclear  Material  (SNM)  identification.  In  order  to  fully 
understand  the  impact  of  various  physics  processes,  as  well  as  to  determine  the  effect  of 
systematic  uncertainties  (such  as  detector  locations),  a  Monte  Carlo  simulation  must  be 
available  that  can  reproduce  the  observed  data.  Here  we  present  a  Monte  Carlo  simula¬ 
tion  that  can  reproduce  characteristic  measurements  of  the  detectors,  such  as  the  charge 
deposited  in  each  scintillator  bar,  as  well  as  the  number  of  bars  hit  per  triggered  event.  The 
MC  results  can  also  be  used  to  create  2D  images  of  cargo,  and  achieve  excellent  agreement 
when  compared  to  data. 
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1  Introduction 


The  Cosmic  Ray  Inspection  and  Passive  Tomography  (CRIPT)  project  is  an  effort  to  con¬ 
struct  a  novel  cosmic  ray  muon  detector  for  the  purposes  of  identifying  Special  Nuclear  Ma¬ 
terial  (SNM).  The  CRIPT  collaboration  has  constructed  a  large-scale  prototype  at  AECL’s 
Chalk  River  Laboratory  to  investigate  the  merits  of  muon  scattering  tomography,  wherein 
the  three-dimensional  density  of  a  scanned  object  is  inferred  by  the  measured  scattering 
of  through-going  charged  particles.  Extremely  dense  materials,  such  as  SNM,  can  then  be 
identified  within  the  large  three-dimensional  scanning  volume  of  the  CRIPT  detector. 

It  is  important  to  evaluate  the  accuracy  of  the  Monte  Carlo  (MC)  simulation.  The  MC  is 
capable  of  turing  off  various  physics  process  and  adding,  moving,  and  removing  material  on 
or  near  the  experiment.  This  allows  analysts  to  identify  features  in  the  data  and  determine 
their  cause.  The  can  lead  to  corrections  to  the  data,  data  selection  criteria,  and  future 
design  (or  modification)  considerations.  The  Monte  Carlo  can  also  be  used  in  run  planing 
and  testing  data  analysis  algorithms  (E.g.  event  tracking,  imaging),  without  the  need  for 
expensive  (in  time  and/or  money)  data  collection. 

1.1  Detector  overview 

The  CRIPT  detector  is  comprised  of  twelve  2x2  m  planes  of  triangular  scintillator  bars. 
Each  plane  contains  121  individual  scintillator  bars,  which  are  2  m  long  and  have  a  tri¬ 
angular  cross-section  with  a  base  of  3.23  cm  and  height  of  1.65  cm.  The  scintillator  bars 
are  arranged  in  a  close-fitting  pattern  of  alternating  right-side-up  and  up-side-down  bars 
such  that  there  are  no  gaps  across  the  plane  surface.  The  scintillated  light  from  each  in¬ 
dividual  bar  is  collected  by  a  wavelength- shifting  optical  fibre  which  is  coupled  to  a  64- 
channel  photo-multiplying  tube  (PMT)  readout  device  (two  PMTs  are  needed  to  read  out 
each  plane). 

Furthermore,  the  CRIPT  detector  is  divided  into  six  layers,  where  each  layer  contains  two 
planes  arranged  orthogonally.  Two  layers  are  located  above  the  imaging  volume  and  are 
referred  to  as  the  Upper  Tracker  (UT),  two  layers  are  located  directly  below  the  imaging 
volume  and  are  referred  to  as  the  Lower  Tracker  (LT),  and  the  final  two  layers  are  located 
below  the  LT  and  are  referred  to  as  the  Spectrometer  (SPEC).  Note  that  this  report  is  mainly 
focused  on  the  optimization  of  the  performance  of  the  UT  and  LT  layers,  as  these  layers  are 
responsible  for  the  tracking  used  for  reconstructive  imaging.  The  SPEC  layers  are  largely 
ignored,  as  the  multiple-scattering  that  occurs  due  to  the  steel  plates  (needed  to  perform  the 
momentum  estimate)  inevitably  degrades  the  overall  SPEC  performance. 


1 


2  Monte  Carlo  Overview 


A  Monte  Carlo  simulation  (MC)  is  a  computer  program  that  calculates  the  physics  pro¬ 
cesses  that  occur  within  a  given  geometry,  provided  a  set  of  start  conditions  and  defined 
primary  particle  characteristics.  The  geometry  includes  the  shape  and  size  of  the  physical 
components,  as  well  as  the  material  that  each  component  is  constructed  of.  To  do  this,  we 
link  to  the  Geant4  software  package[l,  2]  which  provides  a  library  of  physics  process  as 
well  as  framework  to  define  a  geometry.  Initially,  the  current  version  was  geant4.9.6.p02, 
but  a  major  update  occurred  to  geant4. 10.00  during  the  course  of  this  study.  This  change  is 
detailed  in  Section  4,  and  all  results  have  been  updated  to  use  geant4. 10.00. 

The  data  acquisition  (DAQ)  for  the  CRIPT  detector[3]  is  processed  through  a  custom  set 
of  procedures,  through  a  program  called  anaTriangle[4].  This  program  outputs  files  for  the 
analysts  using  a  open-source  analysis  package  called  ROOT[5,  6],  version  5.17.  This  is  a 
standard  analysis  package  in  particle/nuclear  physics. 

The  state  of  the  CRIPT  Monte  Carlo  was  uncertain  as  of  October  2013.  Two  version 
existed,  both  were  partial: 

1.  A  toy  MC  had  some,  but  not  all  detector  components,  but  could  produce  output  that 
could  be  read  by  the  data  analysis  tools, 

2.  A  more  complete  geometry,  but  produced  output  that  could  not  be  analyzed  in  the 
same  way  as  the  data. 

The  problems  with  1.  are  that  we  could  not  count  on  the  output  to  be  accurate,  this  was 
useful  as  a  first-order  check  on  some  features,  but  no  in-depth  studies  could  be  completed. 
2.  had  more  accurate  results  but  could  not  be  process  with  anaTriangle,  so  when  differences 
arose  between  data  and  MC  the  analyst  could  not  be  confident  if  they  due  how  the  MC  was 
processed  or  if  the  MC  really  was  different  than  the  data. 

A  new  MC  was  created  to  replace  both  1.  and  2.  above.  The  goals  of  this  new  MC  were 
to:  have  accurate  detector  components,  optimized  software  design,  and  it’s  output  would  be 
able  to  be  processed  through  anaTriangle,  in  the  same  way  as  the  data.  This  new  MC  would 
eliminate  any  question  to  the  cause  of  features  in  the  MC  vs.  data  comparison.  Figure  1 
was  produced  by  the  new  MC  simulation. 
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Figure  1:  The  MC  geometry  of  the  CRIPT  detector.  The  yellow  bars  at  the  top  are  supports 
for  people  walking  the  top,  the  blue  horizontal  plates  are  steel  flooring  for  walking  (top) 
and  for  cargo  (middle).  The  blue  horizontal  beams  are  the  support  structure.  The  yellow 
plates  near  the  bottom  are  the  iron  slabs  of  the  spectrometer  The  red  plates  are  the  panels 
containing  the  scintillator  bars. 
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3  Monte  Carlo  Geometry 


The  geometry  of  the  detector  is  required  by  the  Geant4  package.  The  structure  from  the 
previous  MCs  were  used  as  a  starting  point.  However,  it  was  discovered  that  there  were 
simplifications  used  in  that  MC  geometry  that  caused  concern.  Such  as: 

•  Components  of  the  geometry  were  overlapping, 

•  parts  had  incorrect  design  (solid  instead  of  hollow), 

•  parts  were  in  wrong  location  (eg.  spectrometer  iron  plates). 

To  fix  this,  Computer  Aided  Design  (CAD)  drawings  were  provided  for  each  component 
of  the  CRIPT  detector.  These  dimensions  and  relative  placement  locations  were  used  in 
the  new  MC.  Figures  1  and  2  show  the  new  detector  geometry,  where  Figure  2  shows  a 
cross-section  of  the  new  detector  geometry.  The  hollow  supports  can  now  be  seen  to  be 
included. 

3.1  Plane  Rotations 

An  important  and  necessary  part  of  the  MC  is  making  sure  the  panel  and  bar  numbering 
schemes  matched  the  data.  This  occurs  implicitly  in  data  and  MC  processing,  and  is  hard  to 
ensure  that  the  data  and  MC  match.  To  do  this  two  comparisons  were  done.  1)  the  charge 
distributions  have  a  slope,  the  coordinates  of  the  MC  had  to  be  set  to  match  between  data 
and  MC,  and  2)  items  were  placed  in  the  cargo  during  running,  and  the  locations  of  these 
could  be  used  to  ensure  the  numbering  scheme  was  implemented  correctly. 

This  became  more  complicated  than  expected,  as  for  historical  reasons  the  y-coordinate 
in  the  MC  was  set  to  be  vertical.  The  data  has  the  vertical  direction  defined  as  the  z- 
coordinate.  This  makes  MC  vs  data  comparisons  difficult  as  there  is  an  axis  direction  flip 
(y  in  MC,  z  in  data).  That  is,  we  have  (unintentionally)  created  a  left-handed  coordinate 
system  (relative  to  data)  in  the  MC.  The  correction  is  to  reverse  the  positive-direction  of  the 
y-bars,  this  looks  backward  when  drawing  the  geometry  but  is  consistent  with  the  physical 
detector,  and  so  is  invisible  to  the  user,  and  to  the  definition  of  the  geometry  placement. 

3.2  Plane  Locations 

Physical  measurements  were  taken  at  Chalk  River  to  determine  where  several  of  the  lower 
layers  and  iron  plates  were  physically  located  in  the  vertical  direction.  However,  these 
measurements  were  taken  after  full  assembly,  so  detector  location  within  packaging  was 
difficult  to  accurately  determine.  An  algorithm  to  fit  for  the  positions  of  the  planes  was 
developed  and  tested  on  the  MC,  then  applied  to  the  data.  This  was  largely  done  as  part 
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Figure  2:  Comparison  of  the  old  (left)  and  new  (right)  MC  geometry. 
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of  Richard  Hydomako’s  work[7].  His  results  were  used  along  with  measurements  for  the 
lower  reference  point  as  inputs  to  the  MC  geometry. 

3.3  Bar  Features 

Each  bar  has  outer  measurement  of  height  of  1.65  cm  and  a  base  of  3.23  cm.  However, 
the  entire  bar  was  not  active.  The  bar  was  constructed  with  an  outer  layer  of  TiO?  which 
had  a  varying  thickness.  This  layer  was  not  sensitive  to  energy  deposition  from  muons  (or 
any  other  particles).  It  is  included  as  a  reflective  material  to  aid  light  collection  from  light 
production  in  the  interior  active  material  (polystyrene  scintillator). 

The  thickness  of  this  had  been  assumed  to  be  0.1  mm.  This  assumption  turned  out  to  be 
an  over  simplification  of  the  parameter.  Figure  3  shows  a  photo  taken  of  two  bars.  We  can 
ignore  the  cloudiness  of  a  bar,  as  this  photo  was  taken  for  another  study.  Since  we  know  the 
outer  dimension  we  set  a  scale  and  measure  the  TiCF  thickness  along  the  walls.  It  measures 
to  be  about  0.3  mm. 


Figure  3:  Photo  of  two  bars  with  the  same  dimensions  as  those  used  in  CRIPT. 

A  more  important  feature  that  can  be  identified  from  Figure  3  is  the  collection  of  TiC>2  in 
the  comers  of  the  bars.  We  can  see  that  there  is  a  much  thicker  layer  of  Ti02  in  the  comers. 
Figure  4  is  the  bar  geometry  used  in  the  MC,  the  additional  Ti02  in  the  comers  can  be  seen. 
This  additional  TiCF  has  a  direct  impact  on  whether  a  single  bar  or  two  bars  are  triggered 
by  a  muon,  as  it  increases  the  cross  section  of  a  single  bar.  The  size  of  deposit  of  TiCF  in 
the  comers  can  be  used  as  a  MC  tuning  parameter  and  is  discussed  in  Section  6.2. 

3.4  Recommendation 

The  plane  rotation  is  fixed,  and  is  well  known  and  tested.  The  plane  locations  are  known 
by  using  data  to  fit  their  positions.  This  position  fit  has  not  been  done  since  the  reassembly 
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Figure  4:  CRIPT  bars  in  the  MC.  The  green  ovals  are  holes  in  the  scintillator  for  the  fibre 
optic  cable. 


of  CRIPT  at  Chalk  River,  and  will  need  to  be  updated  for  new  data. 

The  additional  TiCb  in  the  comers  and  to  a  lesser  extent  the  thickness  of  the  TiC>2  along 
the  walls  will  vary  from  bar-to-bar.  In  the  MC,  this  parameter  has  been  set  to  be  the  same 
for  each  bar,  so  we  are  averaging  over  all  variations.  Due  to  this  implementation,  this 
parameter  can  be  tuned  to  achieve  better  agreement  with  the  data. 
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4  Geant4  version  update 


At  the  start  of  this  contract  period  (October  2013),  GEANT4.9.6.p02  was  the  most  up  to 
date  version.  On  December  6,  2013  the  Geant4  collaboration  released  Geant4. 10.00. 
This  was  a  major  update  that  changed  many  of  the  default  physics  within  the  MC.  It  also 
added  multi-threading  support,  which  is  very  useful  modern  CPUs  have  multiple  cores  that 
can  now  easily  be  used  to  increase  running  speed. 

The  MC  implementation  of  the  muon  multiple  scattering  is  the  most  important  to  an  exper¬ 
iment  like  CRIPT.  Tests  were  done  to  check  how  large  of  an  effect  the  update  would  have 
on  our  results.  To  test  the  muon  multiple  scattering,  we  used  a  20  x  20  x  20  cm  cube  of 
lead.  Mono-energetic  (E=4  GeV)  p  s  were  generated  10  cm  above  the  centre  of  the  block, 
with  a  direction  of  (0,-l,0),  that  is  straight  down  in  the  MC.  The  set  up  and  angle  are  shown 
in  Figure  5.  The  direction  of  the  muon  exiting  the  block  was  then  recorded. 


Figure  5:  A  cartoon  of  the  sep-up  for  Sections  4  and  5.  The  4  GeV  muons  come  from 
the  top  (0,-l,0),  the  red  block  is  a  20  x  20  x  20  cube  of  lead,  and  the  scattering  angle,  0  is 
measured  at  the  bottom. 

This  exiting  direction  was  then  projected  onto  the  xy-plane  (where  y  is  the  vertical  di¬ 
rection).  The  resulting  angle  is  referred  to  as  the  scattering  angle,  0,  and  is  shown  in 
Figure  6.  We  can  see  that  the  width  in  version  GEANT4.9.4  differs  from  GEANT4.9.6  and 
Geant4.  10,  but  the  latter  two  agree.  In  any  case,  the  difference  is  smaller  than  our  detector 
angular  resolution,  so  this  is  not  an  issue  for  our  analysis. 
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Geant4.10 


Figure  6:  A  comparison  of  the  scattering  angle  (xy -projection)  through  20  cm  (cube)  of 
Pb,  for  three  versions  of  Geant4. 

4.1  Recommendation 

The  change  to  the  Geant4  default  settings  and  included  physics  processes  have  not  caused 
a  significant  change  to  the  muon  multiple  scattering  results.  Any  differences  that  we  have 
observed  have  been  significantly  smaller  than  our  detector’s  angular  resolution,  and  there¬ 
fore  can  be  neglected. 
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5  Monte  Carlo  Physics  Lists 


The  Geant4  software  package  requires  that  the  user  define  the  physics  processes  that 
he/she  would  like  to  include  in  the  simulation.  We  set  our  simulation  to  include  all  electro¬ 
magnetic  (EM)  processes,  as  those  are  the  main  ways  in  which  muon  deposit  energy.  We 
should  note  that  we  do  not  include  any  neutron  physics,  as  we  have  assumed  that  neutrons 
cannot  interact  in  a  way  that  would  pass  the  CRIPT  DAQ  trigger.  Including  neutron  physics 
can  cause  a  significant  loss  performance,  as  neutrons  have  long  mean  free  paths. 

The  only  remaining  question  that  we  had  was  if  we  needed  to  include  muon-nuclear  inter¬ 
actions.  To  test  this  we  used  the  same  set-up  as  in  Section  4,  and  turned  the  muon-nuclear 
physics  on/off  in  Geant4.  10.00.  The  scattering  angle  was  measured  with  respect  to  the 
primary  generate  muons  (0,-l,0).  and  the  direction  cosine  was  determined.  The  result  is 
presented  in  Figure  7. 


Figure  7:  A  comparison  of  the  scattering  angle  (direction  cosine)  through  20  cm  (cube)  of 
Pb,  using  GEANT4.1 0.00.  In  one  case  we  include  muon-nuclear  interactions,  in  the  other 
case  we  do  not  include  muon-nuclear  interactions. 

From  Figure  7  we  can  see  no  significant  difference  between  including  or  omitting  the 
muon-nuclear  physics. 

5.1  Recommendation 

As  Figure  7  show  no  significant  difference,  we  choose  to  include  the  muon-nuclear  term, 
as  it  does  not  affect  the  run  performance.  It  may  be  worth  time  looking  at  the  inclusion  of 
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neutron  physics,  if  the  CRIPT  collaboration  is  interested  in  looking  at  neutron  producing 
cargo  or  the  CRIPT  trigger  conditions  set  to  be  more  inclusive  (higher  acceptance). 
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6  Monte  Carlo  vs.  Data 


The  goal  of  any  physics  Monte  Carlo  simulation  is  to  accurately  reproduce  the  physics  in 
the  provided  geometry.  In  the  case  of  CRIPT,  we  have  taken  a  lot  of  data  and  have  many 
figure-of-merit  distributions  that  we  can  produce  in  both  the  data  and  MC.  The  following 
sub-sections  describe  these  figure-of-merit  distributions,  show  some  preliminary  results 
for  the  intrinsic  resolution,  and  finally  present  a  comparison  of  data  and  MC  for  materials 
that  have  been  imaged  in  the  cargo  bay.  Sections  6.5  and  6.6  provide  details  on  how  the 
images  were  made  and  how  we  can  use  the  images  to  identify  the  material  in  the  image, 
respectively. 

6.1  Charge  Collection  and  Distribution 

When  a  muon  interacts  with  the  plastic  scintillator  bar  it  deposits  some  of  its  energy  into 
the  plastic  (^f ).  This  energy  is  converted  to  light,  which  can  travel  and  scatter  within  the 
plastic  and  enter  a  fibre  optic  cable  that  guides  the  photons  to  a  PMT  that  coverts  this 
light  to  a  charge.  This  charge  is  then  proportional  to  the  energy  deposited  by  the  muon. 
The  amount  of  charge  collected  per  unit  energy  deposited  varies  from  bar-to-bar  due  to 
collection  efficiencies  in  the  bar  and  PMT,  as  well  as  light  loss  in  the  bar  and  the  fibre  optic 
cable  (attenuation).  These  bar-to-bar  variations  can  be  calibrated  out  so  that  we  see  the 
same  charge  per  unit  energy  in  all  bars. 

The  mean  charge  in  a  given  bar,  after  calibration,  should  then  be  the  same  for  each  bar  on 
a  given  panel.  This  can  be  seen  for  the  data  in  Figure  8.  The  MC  directly  measures  this 
energy  deposited  and  a  separate  simulation  in  anaTriangle  considers  all  the  imperfections 
in  the  light  collection  and  propagation,  and  electronic  effects  in  the  data.  These  are  also 
calibrated,  and  these  results  are  shown  in  Figure  9.  We  obtain  a  good  agreement  between 
data  and  MC  for  the  mean  charge  distributions. 
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Figure  8:  Charge  distribution  in  the  data  for  a  run  with  an  empty  cargo  bay. 
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Figure  9:  Charge  distribution  in  the  MC  for  a  run  with  an  empty  cargo  bay 
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6.2  Ratio  of  Isolated  Vs.  Adjacent  Bars 

When  a  muon  passes  through  a  detector  panel  it  can  interact  with  one  or  two  (rarely  more) 
triangular  scintillator  bars.  This  ratio  is  effectively  the  cross-section  of  a  given  bar(s).  The 
cross-section  of  a  single  bar  can  be  increase  by  enlarging  the  apex  of  the  triangle  by  adding 
more  TiC>2  (shown  in  Section  3.3)  -  The  TiCF  does  not  produce  photons  when  a  muon 
interacts  within  it.  The  thickness  of  this  TiCb  can  be  measured  from  samples  of  the  bars, 
and  tuned  within  the  measurement  uncertainty. 


The  ratio  of  isolated  to  adjacent  bar  triggers  in  the  data  is  shown  in  Figure  10.  If  we  assume 
in  the  MC  a  thin  layer  of  TiCF,  the  same  as  on  the  flat  sides  of  the  bar,  we  get  ratio  much 
larger  than  the  data,  as  shown  in  Figure  11.  However,  if  we  then  use  the  measured  TiC>2 
thickness  in  the  MC  we  get  a  much  closer  ratio,  as  we  can  see  in  Figure  12. 
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Figure  10:  Using  the  data,  this  is  the  number  of  counts  per  bar  when  either  only  one  bar  is 
triggered  (isolated)  or  two  neighbouring  bars  are  triggered  (adjacent)  for  each  plane. 
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Figure  11:  Using  the  MC  before  accounting  for  extra  TiO 2  in  the  corners  of  the  bar,  this 
is  the  number  of  counts  per  bar  when  either  only  one  bar  is  triggered  (isolated)  or  two 
neighbouring  bars  are  triggered  (adjacent)  for  each  plane. 
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Figure  12:  Using  the  MC  after  accounting  for  extra  T1O2  in  the  corners  of  the  bar,  this 
is  the  number  of  counts  per  bar  when  either  only  one  bar  is  triggered  (isolated)  or  two 
neighbouring  bars  are  triggered  (adjacent)  for  each  plane. 
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6.3  Resolution 


The  intrinsic  resolution  is  the  resolution  of  the  detector  when  you  remove  the  effect  of 
multiple  scattering.  We  believe  the  MC  accounts  for  the  multiple  scattering  properly,  based 
on  the  excellent  agreement  between  MC  and  data  in  all  the  figure  of  merits  that  tested.  This 
means  that  we  define  the  intrinsic  resolution,  Rm,  by  Equation  1: 

R\n  —  -rirue  Vpos  (1) 

where  ,i'tme  is  the  true  location  of  the  interaction  and  ,v'pos  is  the  detected  location  of  the 
interaction.  In  the  case  of  a  single  bar,  this  is  just  the  position  of  that  bar.  If  multiple 
bars  are  triggered  then  the  position  is  charged  weighted  between  the  two  bars.  We  define 
the  true  position,  jctrUe,  by  drawing  a  line  between  the  first  scatter  in  the  panel  and  the  last 
scatter.  We  then  determine  the  coordinate  at  mid-point  of  the  plane  (there  are  12  planes  in 
CRIPT),  and  call  this  the  true  location. 

Figure  13  (top)  shows  the  results  for  Equation  1,  for  the  case  of  isolated  bars  per  trigger  and 
adjacent  bars.  We  can  see  that  we  get  better  resolution  when  we  use  events  that  hit  more 
than  one  bar,  this  makes  sense  as  we  can  get  better  position  information.  Figure  13  (Bot¬ 
tom)  shows  the  same  information  for  resolution,  but  includes  the  location  along  the  length 
of  the  panel.  The  values  in  Figure  13  (Bottom)  are  the  mean  ±  RMS  of  the  distributions. 
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Panel  0  (Corr  =  -0.00547484)  Panel  1  (Coit  =  -0.00428564)  Panel  2  (Corr  =  -0.00489966)  Panel  3  (Corr  =  -0.012214) 


Panel  8  (Corr  =  0.000347822) 


Panel  10  (Corr  =  -0.00379916) 


Panel  11  (Corr  =  0.00370513) 


Figure  13:  Top:  MC  resolution  along  the  length  of  the  bar  for  each  panel.  Bottom:  Pro¬ 
jection  (y)  of  the  top  figure,  divided  between  isolated  bars  and  adjacent  bars. 
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6.4  Cargo 

The  CRIPT  detector  is  designed  to  image  material  that  is  located  in  it’s  cargo  bay.  That  is, 
the  space  located  between  the  UT  and  LT.  Various  items  have  been  imaged,  and  a  selection 
of  these  have  been  added  to  the  MC.  One  item  is  a  depleted  uranium  (DU)  sphere,  with  a 
hole  part  way  through  the  sphere.  This  hole  is  used  to  test  the  resolution  of  the  imaging 
algorithms.  The  remaining  items  that  have  been  included  in  the  MC  are  lOx  lOx  10  cm 
cubes  of  W,  Al,  Cu,  Fe,  Pb,  or  brass.  There  were  multiple  runs  were  different  combinations 
of  these  cubes  were  included  (sets  of  three  or  four).  To  further  test  the  sensitivity  of  our 
imaging  algorithms  various  arrangements  of  lead  bricks  were  used.  Such  as,  arranging 
the  bricks  at  varying  distances  from  each  other  to  test  if  the  resolution  is  good  enough 
to  distinguish  the  individual  bricks  and  bricks  making  up  larger  cubes,  some  with  bricks 
removed  to  test  the  imaging  resolution.  A  55  Gal.  drum  was  filled  with  wax  and  had  a 
W  rod  placed  within  the  wax.  This  would  test  if  we  can  see  a  heavy  material  embedded 
in  a  matrix.  The  bricks  were  also  arranged  to  spell  CRIPT,  this  is  ideal  for  testing  for 
rotations  between  the  data  and  the  MC  through  the  imaging  routines  -  it  also  is  nice  for 
presentations. 

Table  1  lists  the  run  numbers  and  cargo  type  for  selected  runs.  The  comparison  of  the  data 
and  MC  for  each  of  the  runs  are  shown  in  Figures  14  to  20.  The  XY-plane  is  the  horizontal 
plane,  and  the  Z-axis  is  the  projection  of  the  median  scattering  angle  from  voxels  onto  the 
2D  XY-plane  (pixels). 


Run  Numbers 

Cargo 

Data  Triggers 

MC  Triggers 

271 

Pb  (spelling  CRIPT) 

2,164,543 

2,350,364 

277 

W,  Cu,  Al,  brass 

2,035,677 

2,848,998 

302 

W,  Cu,  Fe,  Pb 

8,551,823 

8,403,324 

342 

Phantoms  (Pb  bricks  placed  various  distances  apart) 

4,339,961 

5,684,579 

366 

Two  20  x  20  cm  lead  cubes,  one  is  hollow 

12,218,289 

11,374,846 

482 

55  Gal.  Drum  with  wax  and  W  rod 

11,446,474 

10,346,546 

486 

DU  sphere  (horizontal  cutout) 

16,151,410 

760,540 

N/A 

Fuel  Bundle  (36  steel  rods) 

N/A 

6,659,747 

Table  1:  Run  numbers,  cargo,  and  number  of  data  and  MC  triggers  used  for  images. 

In  Figure  14  we  see  that  we  get  a  very  close  reproduction  of  the  Pb  CRIPT  letters,  although 
the  spacing  seems  a  little  different  in  the  MC  vs.  the  data.  This  is  due  to  the  placement  of 
the  Pb  bricks  in  the  MC,  not  an  imaging  artifact. 

Figure  15  shows  the  arrangement  for  run  277.  We  can  see  that  for  the  same  number  of  trig¬ 
gers  in  the  MC  and  data  the  images  look  very  similar.  A  distinct  feature  in  this  arrangement 
is  that  CRIPT  does  not  have  much  sensitivity  to  Al,  which  is  evident  in  both  the  data  and 
the  MC. 
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Figure  14:  Imaging  run  271.  Data  on  the  left,  MC  on  the  right.  XY  are  the  horizontal 
planes  in  this  image,  Z  is  proportional  to  the  scattering  angle.  Pb  bricks  spelling  CRIPT. 
The  spacing  between  letters  looks  different  between  data  and  MC. 
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Figure  15:  Imaging  run  277.  Data  on  the  left,  MC  on  the  right.  XY  are  the  horizontal 
planes  in  this  image,  Z  is  proportional  to  the  scattering  angle.  From  top  left  to  bottom 
right:  W,  Cu,  Al,  and  brass. 


The  images  in  Figure  16  show  the  arrangement  for  run  302  with  similar  results  to  Figure  15, 
that  is,  that  the  imaging  routine  give  consistent  results  between  data  and  MC. 

For  run  342  we  used  Pb  bricks  with  various  spacing,  the  results  are  shown  in  Figure  17. 
The  goal  is  to  determine  the  ability  of  CRIPT  to  separate  individual  Pb  bricks.  The  spacing 
between  the  bricks  ranges  from  1-12  cm.  In  both  the  MC  and  the  data,  the  smallest  spacing 
distinguishable  is  6  cm.  This  agreement  between  data  and  MC  provides  further  confidence 
that  the  MC  reproduces  the  data  well. 

In  Figure  18,  run  366,  we  test  CRIPTs  sensitivity  to  hollow  section  in  a  Pb  block.  Both 
the  data  and  MC  are  able  to  observe  this  hollow  space  in  their  images.  Like  Figure  17,  this 
provides  a  strong  statement  that  the  MC  agrees  with  the  data. 

Figure  19,  run  482,  shows  the  image  of  a  55  Gal.  drum  with  a  1  x  1  x  8-inch  W  rod  near 


20 


X  (cm)  X  (cm) 


Figure  16:  Imaging  run  302.  Data  on  the  left,  MC  on  the  right.  XY  are  the  horizontal 
planes  in  this  image,  Z  is  proportional  to  the  scattering  angle.  From  top  left  to  bottom 
right:  W,  Cu,  Fe,  and  Pb. 
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Figure  17:  Imaging  run  342.  Data  on  the  left,  MC  on  the  right.  XY  are  the  horizontal 
planes  in  this  image,  Z  is  proportional  to  the  scattering  angle.  Pb  bricks  placed  various 
distances  apart  -  which  gaps  we  are  able  to  distinguish  is  a  measure  of  resolution/precision. 
There  is  an  x-llip  between  MC  and  data. 


the  centre  of  the  drum.  In  the  data  (left)  we  can  start  to  see  the  W  rod,  however  in  the  MC 
(right)  it  is  hard  to  locate  the  rod.  This  is  probably  due  to  the  distribution  of  the  wax  in  the 
data.  The  wax  is  made  up  of  many  square/cylinder  wax  blocks,  resulting  in  only  60-70% 
of  the  volume  being  filled.  The  MC  fills  the  volume  completely  with  paraffin  wax.  This 
difference  doesn’t  seem  to  matter  in  the  imaging,  which  makes  sense,  as  we  saw  about  our 
imaging  is  not  sensitive  enough  to  see  Al,  so  wax  should  not  influence  the  image. 

Finally,  Figure  20  looks  at  a  depleted  uranium  (DU)  sphere,  with  a  hole  bored  into  one 
side.  Both  the  data  and  the  MC  have  trouble  imagine  the  hole.  Even  though  this  may  seem 
disappointing,  the  goal  of  this  study  is  to  produce  a  MC  that  matches  well  to  the  data,  so 
matching  it’s  limitations  evidence  of  a  successful  MC. 
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Figure  18:  Imaging  run  366.  Data  on  the  left,  MC  on  the  right.  XY  are  the  horizontal 
planes  in  this  image,  Z  is  proportional  to  the  scattering  angle.  Two  20  x  20  x  20  cm  lead 
cubes,  one  is  hollow  (x~  80  cm),  the  second  is  solid  (x&  120  cm).  Right  hollow,  left  solid. 
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Figure  19:  Imaging  run  482.  Data  on  the  left,  MC  on  the  right.  XY  are  the  horizontal 
planes  in  this  image,  Z  is  proportional  to  the  scattering  angle.  55  Gal.  drum  filled  with 
wax,  with  a  1  x  1  x  8-inch  W  rod  near  the  centre  of  the  drum. 
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Figure  20:  Imaging  run  486.  Data  on  the  left,  MC  on  the  right.  XY  are  the  horizontal 
planes  in  this  image,  Z  is  proportional  to  the  scattering  angle.  Depleted  uranium  (DU) 
sphere,  with  cutout  in  the  horizontal  direction  (parallel  to  XY-plane). 
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6.4.1  Future  Runs 


Now  that  the  CRIPT  detector  is  located  at  AECL’s  Chalk  River  Laboratories,  it  is  now 
possible  to  image  CANDU  reactor  fuel  bundles.  An  example  bundle  created  with  (solid) 
steel  rods  in  the  same  geometry  as  the  (hollow)  fuel  rods  is  available.  Figure  21  shows  the 
fuel  bundle  geometry  in  the  MC.  Figure  22  show  the  image  of  the  steel  fuel  bundle  for  a 
24  hour  run,  with  the  bundle  in  the  vertical  position. 


After  24  hours  and  zooming  in  on  the  image,  it  is  still  not  possible  to  distinguish  the 
individual  fuel  rods.  Longer  runs  are  currently  underway  at  Chalk  River. 


X  (cm) 

Figure  22:  Imaging  of  steel  MC  fuel  bundle.  0.6  cm 3  bins. 
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6.5  Imaging  Using  Median  vs.  Mean  Scattering  Angle 

When  creating  the  2D  images  in  Section  6.4  we  first  start  with  a  3D  version.  To  create  an 
image  we  use  the  results  from  a  fit  algorithm  known  as  PoCA,  a  simplified  description  of 
how  it  works  is: 


1.  Assume  a  linear  track  from  from  the  upper  tracker 

2.  Assume  a  linear  track  from  the  lower  tracker 

3.  Extrapolate  the  two  lines  from  1.  and  2.  into  the  cargo  area,  until  they  meet. 

4.  The  angle  between  the  lines  from  1.  and  2.  is  the  called  the  scattering  angle. 

5.  The  value  of  this  scattering  angle  is  recorded  in  the  voxel  at  the  intersection  point  in 
the  cargo  bay 

Once  all  of  the  triggered  muons  have  been  processed  each  voxel  can  contain  a  number  of 
entries  (scattering  angles).  Some  examples  of  these  are  shown  in  Figure  23.  Each  voxel  is 
then  represented  by  a  single  number  to  represent  the  distribution.  There  are  two  obvious 
choices:  the  mean  and  the  median.  The  mean  can  be  skewed  by  tails  in  the  distribution, 
however  the  median  is  always  located  with  50%  above  and  below,  so  we  like  to  use  the 
median,  as  it  is  more  robust.  Figure  23  shows  two  lines,  the  mean  and  the  median,  to 
illustrate  this  point. 

To  create  the  2D  images,  we  then  project  (sum  all  bins  along  Z)  the  3D  voxels  onto  the 
XY  plane.  Cuts  can  be  made  to  clean  up  the  image  to  reduce  noise.  The  most  effective  cut 
we  have  found  is  to  check  that  neighbouring  voxels  have  count  rates  above  a  predefined 
threshold,  this  eliminates  isolated  voxels  with  a  high  rate,  a  common  noise  characteristic. 
However,  how  to  apply  this  cut  can  be  tricky. 


There  are  two  way  to  accomplish  this  cut,  first  by  summing  all  the  scattering  angles  (SA) 
in  a  voxel  and  comparing  it  (and  its  neighbours)  to  a  threshold,  as: 

C  =  £SAf,  (2) 

i 


and  second,  by  calculating  the  median  (or  mean)  of  each  voxel,  then  comparing  that  to  a 
threshold: 


SA  ,■ 

Ni 


(3) 


It  turns  out  that  applying  cuts  from  Equations  2  and  3  do  not  give  the  same  result.  This  is 
because  D  has  added  a  function  of  the  number  of  counts  in  the  voxel,  which  varies  from 
voxel-to-voxel  and  thus  creates  a  varying  threshold.  This  creates  more  noise  in  the  image, 
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Figure  23:  Distribution  of  scattering  angles  for  selected  voxels. 


often  presented  as  a  halo  in  the  image.  We  apply  cut  C  in  our  algorithm,  and  then  calculate 
the  median  in  each  voxel.  We  then  project  this  result  to  create  the  2D  image. 

If  we  apply  the  same  cut  (C)  and  compare  the  image  using  the  mean  or  median,  so  that  the 
only  difference  comes  from  this  effect,  we  get  the  images  in  Figure  24.  By  eye,  it  looks 
like  the  figure  on  the  right  (median)  has  a  better  signal-to-noise  ratio  then  the  figure  on  the 
left  (mean).  To  confirm  this  a  ROC  curve  would  have  to  be  created. 
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Figure  24:  Mean  vs.  Median  image  for  run  482.  Mean  on  the  left,  Median  on  the  right. 
XY  are  the  horizontal  planes  in  this  image. 
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6.6  Material  (Cargo)  Identification 

Using  the  cargo  from  Section  6.4,  it  is  possible  to  investigate  the  scattering  angle  distri¬ 
butions  of  various  materials.  Figure  25  shows  the  distributions  for  the  objects  in  run  277 
(Figure  15;  from  top  left  to  bottom  right:  W,  Cu,  Al,  and  brass).  This  figure  is  made  by 
looking  a  region  around  each  cube  in  Figure  15,  and  uses  the  same  number  of  events  in 
both  the  data  and  MC.  I.E.  The  image  of  the  MC  in  Figure  15  was  recalculated  using  only 
the  first  2,035,677  events  to  match  the  data.  The  value  of  each  pixel  (2D  bin)  in  Figure  15 
is  the  scattering  angle.  Each  of  these  pixel  values  in  the  region  containing  the  respective 
cube  is  placed  in  the  corresponding  histogram  in  Figure  25. 


Figure  25:  Distribution  of  median  scattering  angles  for  objects  in  Figure  15.  Data  on  the 
left,  MC  on  the  right.  XY  are  the  horizontal  planes  in  this  image. 

Figure  25  shows  two  sets  of  distributions:  data  on  the  left  and  MC  on  the  right.  Both  show 
the  same  trend,  the  heavier  materials  have  larger  scattering  angles,  as  we  expect.  Note  that 
the  Counts  on  the  y-axis  come  from  the  number  of  bins  in  the  2D  image.  To  determine  an 
uncertainty  we  would  have  to  determine  the  uncertainty  in  the  scattering  angle  distribution 
(such  as  the  ones  shown  in  Figure  23). 

6.7  Summary 

Sections  6.1  and  6.2  provide  us  with  characteristic  distributions  to  aid  in  MC  tuning.  Tun¬ 
ing  the  MC  is  a  necessary  process  to  make  the  MC  replicate  the  data.  The  tuning  provides 
a  way  for  the  MC  to  account  for  parameters  that  are  not  simulated,  such  as  photon  at¬ 
tenuation  lengths  and  DAQ  limitations,  among  others.  The  tuning  should  not  account  for 
inaccuracies  in  the  MC  geometry,  the  geometry  needs  to  be  as  correct  as  possible  to  al¬ 
low  tuning  to  obtain  the  greatest  accuracy  and  precision.  The  tuning  is  carried-out  during 
post-processing  using  anaTriangle. 

Using  the  imaging  algorithm  on  the  MC  gives  similar  results  to  the  images  of  the  data.  We 
can  see  from  Figures  14  to  20  that  we  get  a  good  qualitative  agreement.  Further  studies 
will  be  necessary  to  quantify  this  agreement. 
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Section  6.6  provided  us  with  a  first  look  at  how  CRIPT  can  be  used  to  identify  materials 
in  the  cargo  bay.  This  is  promising,  as  this  is  the  goal  of  CRIPT.  It  is  also  encouraging 
that  the  distributions  make  physical  sense,  that  is  the  heavier  materials  to  correlate  to  larger 
scattering  angles. 
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7  Material  Scattering  and  Detection  due  to 
Muon  Momentum 


Section  6  presented  the  data  vs.  MC  results,  and  how  we  can  use  CRIPT  to  identify  ma¬ 
terials  in  the  cargo  bay.  Here,  I  studied  run  277,  building  on  Section  6.6,  we  can  examine 
how  the  muon  momentum  affects  the  scattering  in  the  cargo  material.  Figure  26  shows  the 
muon  momentum  distribution  for  events  that  cause  a  trigger  in  CRIPT. 


Figure  26:  Distribution  of  triggered  muons  for  run  211. 

Figures  15  and  25  show  the  2D  image  reconstruction  and  material  scattering  distribution 
for  run  277 .  The  following  figures  will  show  subsections  of  these,  cutting  on  the  gener¬ 
ated  momentum  of  the  muon,  p jU.  I  should  note  here,  that  this  can  only  be  done  for  the 
MC,  as  it  allows  us  to  know  the  exact  generated  muon  momentum.  Figure  27  shows  the 
generated  muon  momentum  for  the  range  0<  <500  MeV,  Figure  28  shows  the  gener¬ 

ated  muon  momentum  for  the  range  500<  pu  <1000  MeV,  Figure  29  shows  the  generated 
muon  momentum  for  the  range  1000<  pu  <2000  MeV,  Figure  30  shows  the  generated 
muon  momentum  for  the  range  2000 <  p ^  <3000  MeV. 

We  can  see  that  for  500  MeV<  pLI  that  we  get  no  information.  This  is  mostly  likely  due 
to  an  excess  of  scattering  and  absorption  of  muons  that  interact  in  the  material,  due  to  the 
low  energy  of  the  muon.  At  higher  energy,  those  with  pu  >2000  MeV  we  see  also  see 
obtain  no  information  about  the  material.  This  is  due  to  the  high  energy  muons  scattering 
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very  little  in  the  material.  This  study  informs  us  that  muons  with  a  incoming  momentum 
of  500 <  pjj  <2000  MeV  are  the  important  events  from  CRIPT. 


0<p<500  MeV 


X  (cm)  (Projection)  Median  Scattering  (PoCA)  Angle  [mRad] 


Figure  27:  Distribution  of  median  scattering  angles  for  objects  in  Run  211  for  the  momen¬ 
tum  range:  0<  p M  <500  MeV.  Image  on  the  left,  median  scattering  angle  on  the  right.  XY 
are  the  horizontal  planes  in  this  image. 


500<p<1000  MeV 


X  (cm)  (Projection)  Median  Scattering  (PoCA)  Angle  [mRad] 


Figure  28:  Distribution  of  median  scattering  angles  for  objects  in  Run  211  for  the  momen¬ 
tum  range:  500<  <  /  000  MeV.  Image  on  the  left,  median  scattering  angle  on  the  right. 

XY  are  the  horizontal  planes  in  this  image. 
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1000<p<2000  MeV 


Figure  29:  Distribution  of  median  scattering  angles  for  objects  in  Run  211  for  the  momen¬ 
tum  range:  1000<  <2000  MeV.  Image  on  the  left,  median  scattering  angle  on  the  right. 

XY  are  the  horizontal  planes  in  this  image. 


p  2000-3000  GeV 


x  (cm)  (Projection)  Median  Scattering  (PoCA)  Angle  [mRad] 


Figure  30:  Distribution  of  median  scattering  angles  for  objects  in  Run  211  for  the  momen¬ 
tum  range:  2000<  pLI  <3000  MeV.  Image  on  the  left,  median  scattering  angle  on  the  right. 
XY  are  the  horizontal  planes  in  this  image. 
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8  Conclusion 


This  work  has  created  a  Monte  Carlo  simulation  that  contains  all  the  relevant  physics  pro¬ 
cess  necessary  for  muon  scattering,  and  includes  all  the  detector  and  structure  geometry  as 
provided  by  the  CAD  drawing  used  for  construction.  The  MC  can  produce  output  that  can 
be  read  and  processed  by  the  same  software  tools  (anaTriangle)  that  process  the  data  from 
the  DAQ  system. 

With  proper  tuning  we  achieve  excellent  agreement  between  the  MC  and  data  for  character¬ 
istic  distributions.  This  includes  the  reproduction  of  the  scattering  angle  of  objects  placed 
in  the  cargo  bay.  The  ability  of  the  MC  to  reproduce  the  data  allows  for  further  testing  our 
imaging  algorithms  and  sensitivity  to  various  material,  such  as  SNM,  without  spending  the 
time  and  money  collecting  data.  The  MC  is  an  excellent  tool  to  use  to  determine  which 
samples  we  should  examine,  and  it  aids  scheduling  and  run  planing. 

The  MC  can  also  be  used  for  planning  future  experiments,  such  as  expanding  CRIPT  to 
include  side  detectors  to  accept  more  horizontal  muons  and  obtain  larger  sensitivity  for  3D 
images.  It  can  also  be  used  to  calculate  the  sensitivity  of  scaled  down  muon-tomography 
experiments. 
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